##################################################################
# /* Author: Gautam Nair and Nicholas Sambanis */
# /* Violence Exposure and Ethnic Identification: Evidence from Kashmir */
# /* Table 15: Heterogeneous Effects of Violence Exposure on Principal Components Index of National Identification*/


##################################################################
# loading packages
##################################################################
library(Hmisc)
library(foreign)
library(sandwich)
library(lmtest)
library(numDeriv)
library(stargazer)
library(ggplot2)
library(plyr)
library(gridExtra)
library(dplyr)
library(scales)

##################################################################

rm(list=ls())

library(stargazer)

data.working <- read.csv("nei_kashmir_replication_dataset.csv")

ind.var <- c("group2_violence")

data.working <- data.working[data.working$group1_control==1 | data.working$group2_violence==1,]

vars <- c("median_age",               
				"gender_female",		
				"education_class12",
				"religion_sunni",
				"main_earner",
				"high_index_income_wealth1",
				"occupation_govt",
				"pc_high_national_id_acc"
)

treatvar <- c("group2_violence")

data.working$i_median_age <- data.working$group2_violence*data.working$median_age
data.working$i_gender_female <- data.working$group2_violence*data.working$gender_female
data.working$i_education_class12 <- data.working$group2_violence*data.working$education_class12
data.working$i_religion_sunni <- data.working$group2_violence*data.working$religion_sunni
data.working$i_main_earner <- data.working$group2_violence*data.working$main_earner
data.working$i_high_index_income_wealth1 <- data.working$group2_violence*data.working$high_index_income_wealth1
data.working$i_occupation_govt <- data.working$group2_violence*data.working$occupation_govt
data.working$i_pc_high_national_id_acc <- data.working$group2_violence*data.working$pc_high_national_id_acc

ivars <- c("i_median_age",               
	"i_gender_female",		
	"i_education_class12",
	"i_religion_sunni",
	"i_main_earner",
	"i_high_index_income_wealth1",
	"i_occupation_govt",
	"i_pc_high_national_id_acc"
)

spec1 <- c(treatvar)
spec2 <- c(treatvar, vars)
spec3 <- c(treatvar, vars, ivars)

my.formula <-paste("pc_index_id_no_ind_pak",'~', paste(spec1, collapse= ' + '))
temp.model.1 <- lm(my.formula, data=data.working)
temp.se.1 <- coeftest(temp.model.1, vcov=vcovHC(temp.model.1,type="HC2"))[,2]

my.formula <-paste("pc_index_id_no_ind_pak",'~', paste(spec2, collapse= ' + '))
temp.model.2 <- lm(my.formula, data=data.working)
temp.se.2 <- coeftest(temp.model.2, vcov=vcovHC(temp.model.2,type="HC2"))[,2]

my.formula <-paste("pc_index_id_no_ind_pak",'~', paste(spec3, collapse= ' + '))
temp.model.3 <- lm(my.formula, data=data.working)
temp.se.3 <- coeftest(temp.model.3, vcov=vcovHC(temp.model.3,type="HC2"))[,2]

varlabels <- c(
"Violence Exposure", 
"Age Above Median", 
"Female",
"Completed High School", 
"Sunni Muslim", 
"Main Earner", 
"Income and Wealth Above Median", 
"Government Occupation",
"Constituency National Identification Above Median", 
"Age x Violence",
"Female x Violence",
"High School x Violence",
"Sunni Muslim x Violence",
"Main Earner x Violence",
"Income and Wealth x Violence", 
"Government Occupation x Violence",
"Constituency National Identification x Violence",
"Constant")
spectitle <- c("Heterogeneous Treatment Effects of Violence Exposure on Identification Index (Principal Components)")
output.file <- c("tf_t_16_cate_violence_pc_index")
temptype <- c("latex", "text")
tempext <- c(".tex", ".txt")

for(q in 1:length(temptype)){
	temp.output <- paste(output.file, tempext[q], sep="")

 stargazer(temp.model.1, temp.model.2, temp.model.3, se=list(temp.se.1, temp.se.2, temp.se.3),
         title= spectitle,
          out=temp.output,
          no.space=TRUE, 
          model.numbers= TRUE,
         notes.append=FALSE,
          align=FALSE,
         covariate.labels=varlabels,
          header= FALSE,
          font.size="small",
          model.names=FALSE,
          digits=2,
          star.char=c("**", "***"),
          star.cutoffs=c(0.05, 0.01),
          omit.stat = c("rsq", "f","adj.rsq"),
          column.labels  = NULL,
          dep.var.caption  = "",
          dep.var.labels.include=FALSE,
         df=FALSE,
         omit.table.layout="n",
         nobs=TRUE,
         type= temptype[q]
        # order= c(spec3, "(Intercept)")
)
}
